function x_new = rk4_lorenz(x0, dt, sigma, rho, beta)
    % rk4_lorenz - 使用四阶 Runge-Kutta 方法对洛伦兹系统进行单步更新
    %
    % 输入:
    %   y0    - 初始状态向量 [u; v; w]（列向量或可以转换为列向量）
    %   dt    - 更新步长
    %   sigma - 洛伦兹系统参数 sigma
    %   rho   - 洛伦兹系统参数 rho
    %   beta  - 洛伦兹系统参数 beta
    %
    % 输出:
    %   y_new - 单步更新后的状态向量 [u; v; w]
    
    % 确保 y0 为列向量
    x0 = x0(:);
    
    k1 = lorenz(x0, sigma, rho, beta);                         % 第一阶段（k1）
    k2 = lorenz(x0 + 0.5*k1, sigma, rho, beta);                % 第二阶段（k2）
    k3 = lorenz(x0 + 0.5*k2, sigma, rho, beta);                % 第三阶段（k3）
    k4 = lorenz(x0 + k3, sigma, rho, beta);                    % 第四阶段（k4）

    % 更新状态
    x_new = x0 + dt/6*(k1 + 2*k2 + 2*k3 + k4);       % 使用RK4公式更新状态
end


function dXdt = lorenz(X, sigma, rho, beta)
    % lorenz - 计算洛伦兹系统状态导数
    %
    % 输入:
    %   X     - 状态向量 [u; v; w]
    %   sigma, rho, beta - 洛伦兹系统参数
    %
    % 输出:
    %   dXdt  - 状态导数向量 [du; dv; dw]

    % 解构状态变量
    u = X(1);
    v = X(2);
    w = X(3);
    
    %根据论文方程（11）计算洛伦兹系统的导数
    du = sigma * (v - u);
    dv = rho * u - v - u * w;
    dw = u * v - beta * w;
    
    dXdt = [du; dv; dw];
end
